Part 2 Fire Selection

2.1 Set Up

2.1.1 Libraries

library(tidyverse)
library(terra)
library(sf)
library(mapview)
library(raster)
# library(rgeos)
library(lubridate)
library(ggplot2)
library(exactextractr)
library(gridExtra)
library(knitr)
# library(rasterVis)
library(RColorBrewer)
library(spData)
library(forcats)
library(cowplot)

2.1.2 USDA National Forest Type Group Dataset

Conifer Forest Type Groups: Douglas-Fir, Fir-Spruce-Mountain Hemlock, Lodgepole Pine

# forest type groups and key
conus_forestgroup <- rast('data/forest_type/conus_forestgroup.tif')
forest_codes <- read_csv('data/forest_type/forestgroupcodes.csv')

# set crs
crs = crs(conus_forestgroup)

2.1.3 EPA level-3 Ecoregions

Canadian Rockies, Idaho Batholith, Middle Rockies, Columbian Mountains - Northern Rockies

# level 3 ecoregions
l3eco <- st_read('data/ecoregion/us_eco_l3.shp') %>%
  st_transform(., crs=crs)

# select northern rocky mountains from level3 ecoregions
eco_select <- l3eco %>% 
  filter(NA_L3NAME %in% c('Canadian Rockies','Columbia Mountains/Northern Rockies','Middle Rockies','Idaho Batholith'))

2.1.4 Mapping

2.1.4.1 Ecoregions

# mapview
palette <- brewer.pal(18,"YlGn")
palette[1] <- rgb(255, 255, 255, maxColorValue=255, alpha=1)
mapview(eco_select,na.color=palette[1],legend=TRUE)

2.1.4.2 Forest Type Groups

# convert raster values to factors
forestgroup_eco <- crop(conus_forestgroup,eco_select) %>% 
  mask(.,eco_select) %>% 
  as.factor()

# add a labels for forest type code 
group_levels <- levels(forestgroup_eco)[[1]]
group_levels[["forest_type"]] <- c("0: Unforested","120: Spruce/Fir","180: Pinyon/Juniper","200: Douglas-fir","220: Ponderosa Pine","240: Western White Pine","260: Fir/Spruce/Mountain Hemlock","280: Lodgepole Pine","300: Hemlock/Sitka Spruce","320: Western Larch","360: Other Western Softwood","370: California Mixed Conifer","400: Oak/Pine","500: Oak/Hickory","700: Elm/Ash/Cottonwood","900: Aspen/Birch","920: Western Oak","950: Other Western Hardwoods")

levels(forestgroup_eco) <- group_levels

# mapview
mapview(raster(forestgroup_eco), col.regions=palette,na.color=palette[1],legend=TRUE)

2.2 Define Fire Parameters

2.2.2 Group Adjacent Fires

# function to group adjoining fire polygons to ensure contiguous high-severity patches
group_fires <- function(mtbs_year) {

  # join the polygons with themselves, and remove those that do not join with any besides themselves
  combined<- st_join(mtbs_year, mtbs_year, join=st_is_within_distance, dist = 180, left = TRUE,remove_self = TRUE) %>% 
    drop_na(Event_ID.y)%>% 
    dplyr::select(Event_ID.x,Event_ID.y)
  
  if(nrow(combined)>=1){ # if there are overlaps for this years fires...
    
    # partition data into that that has overlap, and that that does not
    overlap <- mtbs_year %>%
      filter(Event_ID %in% combined$Event_ID.x)
    no_overlap <- mtbs_year %>%
      filter(!(Event_ID %in% combined$Event_ID.x))
    
    print(paste0("there are ",nrow(overlap)," overlapping polygons"))
    
    # join all overlapping features, and buffer to ensure proper grouping
    overlap_union <- st_union(overlap) %>%
      st_buffer(190)
    
    # break apart the joined polygons into their individual groups
    groups <- st_as_sf(st_cast(overlap_union ,to='POLYGON',group_or_split=TRUE)) %>%
      mutate(year = mean(mtbs_year$year),
             Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
      rename(geometry = x)
    
    print(paste0("polygons formed into ",nrow(groups)," groups"))
    
    # join back with original dataset to return to unbuffered geometry
    grouped_overlap <- st_join(overlap,groups,left=TRUE)
    
    # arrange by the new grouping
    joined_overlap_groups <- grouped_overlap %>%
      group_by(Fire_ID) %>%
      tally()%>%
      st_buffer(1) %>%
      dplyr::select(Fire_ID) %>%
      mutate(year = mean(mtbs_year$year))
    
    # add new ID to the freestanding polygons
    no_overlap_groups <- no_overlap %>%
      mutate(Fire_ID = str_c("Fire_",nrow(groups)+c(1:nrow(no_overlap)),"_",year)) %>%
      dplyr::select(Fire_ID,year)
    
    # join the new grouped overlap and the polygons without overlap
    fires_export <- rbind(joined_overlap_groups,no_overlap_groups)
    return(fires_export)
    
    } else { # if there are no overlaps for this year...
      
      print("no overlapping polygons")
      
      fires_export <- mtbs_year %>%
        mutate(Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
        dplyr::select(Fire_ID,year)
      
      return(fires_export)
  }
}
# group adjacent polygons within each fire year
fires_88 <- group_fires(mtbs_select %>%  filter(year == 1988))
## [1] "there are 22 overlapping polygons"
## [1] "polygons formed into 7 groups"
fires_89 <- group_fires(mtbs_select %>%  filter(year == 1989))
## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_90 <- group_fires(mtbs_select %>%  filter(year == 1990))
## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_91 <- group_fires(mtbs_select %>%  filter(year == 1991))
## [1] "no overlapping polygons"
# join each fire year, filter by area
mtbs_grouped <- rbind(fires_88,fires_89,fires_90,fires_91)%>%
  mutate(area_ha = as.numeric(st_area(geometry))/10000,
         area_acres = area_ha*2.471)

2.2.3 Select Fires by Ecoregion and Forest Type

# assign ecoregion and proportions of forest type to each fire polygon
fires_join <- st_join(mtbs_grouped,eco_select,join=st_intersects,left=FALSE,largest=TRUE) %>% 
  left_join(., exact_extract(conus_forestgroup,mtbs_grouped, append_cols = TRUE, max_cells_in_memory = 3e+08, 
                             fun = function(value, coverage_fraction) {
                               data.frame(value = value,
                                          frac = coverage_fraction / sum(coverage_fraction)) %>%
                                 group_by(value) %>%
                                 summarize(freq = sum(frac), .groups = 'drop') %>%
                                 pivot_wider(names_from = 'value',
                                             names_prefix = 'freq_',
                                             values_from = 'freq')}) %>%
              mutate(across(starts_with('freq'), replace_na, 0)))
 
# remove unnecessary columns, cleanup names
# filter to ensure fire polygons are at least 25% type of interest
fires <- fires_join %>% 
  dplyr::select("Fire_ID","year","area_ha","area_acres","US_L3NAME","freq_0","freq_200","freq_220","freq_260","freq_280") %>% 
  rename("ecoregion" = "US_L3NAME",
         "freq_df"="freq_200",
         "freq_pp"="freq_220",
         "freq_fs"="freq_260",
         "freq_lpp"="freq_280") %>% 
  mutate(freq_allother = 1-(freq_0 + freq_df+freq_pp+freq_fs+freq_lpp),
         freq_forested = 1- freq_0,
         freq_ideal = freq_df+freq_fs+freq_lpp)%>% 
  mutate(across(starts_with('freq'), round,2))%>% 
  filter(freq_ideal > 0.25)

2.2.4 Select Fires by Burn Severity

# import all mtbs rasters via a list
rastlist <- list.files(path = "data/mtbs", pattern='.tif', all.files=TRUE, full.names=TRUE)
allrasters <- lapply(rastlist, raster)
names(allrasters) <- str_c("y", str_sub(rastlist,22,25))

# create empty dataframe
severity_list <- list()

# loop through mtbs mosasics for 1988-1991
# extract mtbs burn severity raster for all selected fires
# calculate burn severity percentages for each fire
for (i in names(allrasters)){
  mtbs_year <- allrasters[[i]]
  fire_year <- filter(fires, year==str_sub(i,2,5)) 
  raster_extract <- exact_extract(mtbs_year,fire_year, max_cells_in_memory = 3e+09,coverage_area=TRUE)
  names(raster_extract) <- fire_year$Fire_ID 
  
  output_select <- bind_rows(raster_extract, .id = "Fire_ID")%>%
    group_by(Fire_ID , value) %>%
    summarize(total_area = sum(coverage_area)) %>%
    group_by(Fire_ID) %>%
    mutate(proportion = total_area/sum(total_area))%>% 
    dplyr::select("Fire_ID","value","proportion") %>% 
    spread(.,key="value",value = "proportion")
  
  severity_list[[i]] <- output_select
}

# combine extracted raster datasets
severity_df <- do.call(rbind, severity_list) 

# join burn severity % to fires polygons
# fix naming
# filter dataset for 500 acres high severity
fires_severity <- left_join(fires,severity_df,by="Fire_ID")%>% 
  rename(noburn= "1",lowsev = "2", medsev = "3", highsev = "4",regrowth = "5", error = "6") %>% 
  dplyr::select(- "NaN",-"regrowth",-"error") %>% 
  mutate(highsev_acres = area_acres*highsev)%>% 
  filter(highsev_acres > 500)

2.2.5 Clean Up Dataset

# get the most common forest type within each polygon
fires_select <- fires_severity %>%
  left_join(.,exact_extract(conus_forestgroup,fires_severity, 'mode', append_cols = TRUE, max_cells_in_memory = 3e+08)) 

fires_select$mode <- as.factor(fires_select$mode)

fires_select <- fires_select %>% 
    mutate(fire_foresttype = case_when(mode==200 ~ "Douglas-Fir",
                                       mode==220 ~ "Ponderosa",
                                       mode==260 ~ "Fir-Spruce",
                                       mode==280 ~ "Lodegepole Pine",
                                       TRUE ~ "Other"),
           Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year))
# join the grouped fires back to original mtbs boundaries
fires_mtbs <- st_join(mtbs_select,fires_select,left=FALSE,largest=TRUE) %>% 
  filter(year.x==year.y)%>% 
  dplyr::select("Event_ID","Incid_Name","Fire_ID","Ig_Date","year.y","state","BurnBndAc","ecoregion") %>% 
  rename(year= year.y)

2.3 Final Fire Dataset

2.3.1 Dataset Overview

full_dataset <- fires_select %>% 
  st_drop_geometry() %>% 
  dplyr::select(Fire_ID,year,ecoregion,fire_foresttype,area_acres,highsev) %>% 
  mutate(highsev = round(highsev,2),
         area_acres = round(area_acres,0))  

kable(full_dataset, 
      align = 'c',
      padding = 1,
      col.names = c("Fire ID", "Year", "Ecoregion", "Majority Forest Type","Area (acres)", "High Severity %"),
      caption = "High-Severity Conifer-Dominated Fires 1988-1991")
Table 2.1: High-Severity Conifer-Dominated Fires 1988-1991
Fire ID Year Ecoregion Majority Forest Type Area (acres) High Severity %
Fire_1_1988 1988 Middle Rockies Fir-Spruce 342005 0.27
Fire_2_1988 1988 Middle Rockies Lodegepole Pine 777690 0.28
Fire_3_1988 1988 Middle Rockies Lodegepole Pine 448911 0.21
Fire_4_1988 1988 Idaho Batholith Douglas-Fir 5651 0.16
Fire_5_1988 1988 Idaho Batholith Fir-Spruce 11945 0.23
Fire_6_1988 1988 Idaho Batholith Douglas-Fir 50666 0.07
Fire_7_1988 1988 Middle Rockies Lodegepole Pine 167870 0.50
Fire_8_1988 1988 Idaho Batholith Douglas-Fir 25889 0.02
Fire_9_1988 1988 Idaho Batholith Douglas-Fir 8312 0.17
Fire_10_1988 1988 Canadian Rockies Lodegepole Pine 42492 0.52
Fire_11_1988 1988 Idaho Batholith Fir-Spruce 45075 0.43
Fire_12_1988 1988 Middle Rockies Lodegepole Pine 5633 0.24
Fire_13_1988 1988 Idaho Batholith Douglas-Fir 4962 0.25
Fire_14_1988 1988 Middle Rockies Other 35864 0.46
Fire_15_1988 1988 Idaho Batholith Fir-Spruce 19499 0.29
Fire_16_1988 1988 Idaho Batholith Fir-Spruce 5626 0.09
Fire_17_1988 1988 Idaho Batholith Fir-Spruce 12746 0.40
Fire_18_1988 1988 Middle Rockies Other 13108 0.48
Fire_19_1988 1988 Idaho Batholith Fir-Spruce 7241 0.27
Fire_20_1988 1988 Idaho Batholith Fir-Spruce 6559 0.13
Fire_21_1988 1988 Middle Rockies Fir-Spruce 3113 0.17
Fire_22_1988 1988 Middle Rockies Other 29233 0.16
Fire_23_1988 1988 Middle Rockies Other 1363 0.41
Fire_24_1988 1988 Idaho Batholith Douglas-Fir 48136 0.31
Fire_25_1988 1988 Northern Rockies Douglas-Fir 8089 0.25
Fire_26_1988 1988 Middle Rockies Douglas-Fir 8588 0.56
Fire_27_1988 1988 Northern Rockies Douglas-Fir 11403 0.30
Fire_28_1988 1988 Northern Rockies Douglas-Fir 21854 0.25
Fire_29_1988 1988 Canadian Rockies Douglas-Fir 33844 0.27
Fire_30_1988 1988 Northern Rockies Douglas-Fir 1909 0.31
Fire_31_1988 1988 Middle Rockies Other 6282 0.47
Fire_32_1989 1989 Idaho Batholith Fir-Spruce 13334 0.15
Fire_33_1989 1989 Middle Rockies Lodegepole Pine 3298 0.37
Fire_34_1989 1989 Idaho Batholith Douglas-Fir 4928 0.38
Fire_35_1989 1989 Idaho Batholith Douglas-Fir 47680 0.02
Fire_36_1989 1989 Idaho Batholith Fir-Spruce 2486 0.30
Fire_37_1989 1989 Idaho Batholith Fir-Spruce 5566 0.30
Fire_38_1989 1989 Idaho Batholith Fir-Spruce 7443 0.28
Fire_39_1989 1989 Idaho Batholith Fir-Spruce 6786 0.26
Fire_40_1989 1989 Idaho Batholith Douglas-Fir 8733 0.12
Fire_41_1989 1989 Idaho Batholith Fir-Spruce 1615 0.49
Fire_42_1989 1989 Idaho Batholith Lodegepole Pine 2488 0.33
Fire_43_1989 1989 Idaho Batholith Fir-Spruce 3081 0.27
Fire_44_1990 1990 Middle Rockies Douglas-Fir 2535 0.45
Fire_45_1990 1990 Idaho Batholith Fir-Spruce 3418 0.53
Fire_46_1990 1990 Idaho Batholith Douglas-Fir 2249 0.49
Fire_47_1990 1990 Middle Rockies Lodegepole Pine 2763 0.41
Fire_48_1990 1990 Middle Rockies Other 13461 0.19
Fire_49_1991 1991 Middle Rockies Douglas-Fir 6978 0.31
Fire_50_1991 1991 Middle Rockies Douglas-Fir 3097 0.18
Fire_51_1991 1991 Middle Rockies Fir-Spruce 6995 0.14
Fire_52_1991 1991 Idaho Batholith Douglas-Fir 1186 0.55
Fire_53_1991 1991 Northern Rockies Fir-Spruce 7095 0.39
Fire_54_1991 1991 Northern Rockies Fir-Spruce 2478 0.51

2.4 Mapping

2.4.1 Selected fires by year

# plot
mapview(fires_select, zcol = "year")

2.4.2 Selected fires by majority forest type

# plot
mapview(fires_select, zcol = "fire_foresttype")

2.4.3 Final Fires

fire_names <- c("Fire_1_1988","Fire_2_1988","Fire_3_1988","Fire_4_1988","Fire_7_1988","Fire_9_1988","Fire_10_1988","Fire_11_1988","Fire_12_1988","Fire_13_1988","Fire_14_1988","Fire_15_1988","Fire_16_1988","Fire_18_1988","Fire_19_1988","Fire_20_1988","Fire_22_1988","Fire_23_1988","Fire_25_1988","Fire_26_1988","Fire_28_1988","Fire_29_1988","Fire_31_1988","Fire_32_1989","Fire_33_1989","Fire_35_1989","Fire_38_1989","Fire_41_1989","Fire_42_1989","Fire_48_1990","Fire_49_1991","Fire_50_1991","Fire_51_1991","Fire_54_1991")

2.4.4 Dynamic Map

states49 <- (us_states["NAME"])%>% st_transform(3488)

mapview(raster(forestgroup_eco), col.regions=palette,na.color="transparent",legend=FALSE) +
  mapview(fires_select %>% filter(Fire_ID %in% fire_names),col.regions="black",alpha.regions=100) +
  mapview(st_union(eco_select), col.regions = "black", color = "black", alpha.regions = 0,lwd=2.5) +
  mapview(states49,alpha.regions=0,legend=FALSE, color = "black")
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 16837548 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  16837548 '

2.4.5 Static Map

crs = 4388
borders <- rnaturalearth::ne_states(c("united states of america", "canada")) %>% st_as_sf()
eco_border <- st_union(eco_select) %>% st_transform(crs = crs) 
conus_forestgroup <- rast('data/forest_type/conus_forestgroup.tif')

forestgroup_map <- crop(conus_forestgroup,st_buffer(eco_select,290000))%>% 
  as.factor() %>% 
  aggregate(fact=3,fun = "modal") %>% 
  project(.,eco_border)
## 
|---------|---------|---------|---------|
=========================================
                                          
## Warning: [project,SpatRaster] argument y (the crs) should be a character value
forestgroup_df <- as.data.frame(forestgroup_map,xy = T) %>% 
  mutate(ftype = case_when(label == 0 ~ "Unforested",
                           label == 200 ~ "Douglas-Fir",
                           label == 220 ~ "Ponderosa Pine",
                           label == 280 ~ "Lodgepole Pine",
                           label == 260 ~ "Fir-Spruce",
                           TRUE ~ "Other"))

forestgroup_df$ftype <- fct_relevel(forestgroup_df$ftype , c("Douglas-Fir","Lodgepole Pine","Fir-Spruce","Ponderosa Pine", "Other", "Unforested"))

palette <- brewer.pal(4,"YlGn")

inset <-ggplot() +
  geom_sf(data = states49 %>% st_transform(crs = crs), fill = NA,color=alpha("#525258"),lwd=.4) +
  geom_sf(data = st_union(eco_select) %>% st_transform(crs = crs),color = "black",fill="black",lwd=.4) +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

mainplot <- ggplot() +
  geom_raster(data = forestgroup_df,aes(x,y,fill=as.factor(ftype)),alpha = .8)+
  scale_fill_manual(values = c("#C3E699","#68CE66","#238443","#8BB62F","#D8E895","white"))+ 
  geom_sf(data = borders %>% st_transform(crs = crs), fill = NA,color=alpha("black",0.14),lwd=.05) +
  geom_sf(data = fires_select %>% filter(Fire_ID %in% fire_names) %>% st_transform(crs = crs),fill = "black", color = "black") +
  geom_sf(data = st_union(eco_select) %>% st_transform(crs = crs), fill = NA, color=alpha("black",0.7),lwd=.4,linetype = "F1") +
  scale_x_continuous(limits = c(-119.4952-.5,-103.21640-1.68), expand = c(0, 0)) +
  scale_y_continuous(limits = c(41.89902-.5,49.00152+.5), expand = c(0, 0)) +
  labs(fill = "Forest Type", x = "Longitude", y = "Latitude") +
   theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

ggdraw() +
  draw_plot(mainplot) +
  draw_plot(inset,
    height = 0.2,
    x = .39,
    y = 0.15
  )
## Warning: Removed 2241609 rows containing missing values (`geom_raster()`).

2.5 Export Data

2.5.1 Final Cleanup for Export

# reformat and project
fires_export <- fires_select %>% 
  mutate(year = as.integer(year)) %>% 
  st_transform(., crs="EPSG:4326")

mtbs_export <- fires_mtbs %>% 
  mutate(year = as.integer(year)) %>% 
  st_transform(., crs="EPSG:4326")

2.5.2 Export

# st_write(fires_export, "data/fire_boundaries/", "fires_export.shp", driver = 'ESRI Shapefile')
# st_write(mtbs_export, "data/fire_boundaries/", "mbts_export.shp", driver = 'ESRI Shapefile')